/* [基本尺寸] */
// 叶片数量
blades = 5; // [1,2,3,4,5,6,7,8]

// 螺旋桨直径(毫米)
propdia = 100; // [30:1:300]

// 叶片螺距 - 通常为叶片长度的75%到133%(-1表示自动)
bladepitch = -1; // [-1:1:200]

// 轮毂直径(毫米)
hubdia = 10; // [3:0.5:30]

// 安装孔直径(毫米)
holedia = 2; // [1:0.1:29]

/* [叶片形状] */
// 最大弦长与叶片长度的比例
maxchordfrac = 0.2; // [0.1:0.05:0.9]

// 轮毂长度(毫米)
hublen = 10; // [3:0.5:50]

// 椭圆轮廓长度系数(>1)
elenfrac = 1.05; // [1.001:0.001:2]

// 旋转方向(1=顺时针, -1=逆时针)
dir = -1; // [-1:2:1]

// 中心线位置(0=前缘, 1=后缘)
centerline = 0.3; // [0:0.1:1]

// 叶片扫掠角(度/毫米)
angle_sweep = 0; // [-5:0.1:5]

/* [翼型属性] */
// 后缘厚度(毫米)
te_thickness = 0.4; // [0.2:0.1:2]

// 根部过渡厚度(毫米)
fairing = 3; // [0:0.5:10]

/* [NACA翼型] */
// 根部NACA代码
naca_root = 8430; // [0:9999]

// 中部NACA代码 
naca_mid = 6412; // [0:9999]

// 尖部NACA代码
naca_tip = 3412; // [0:9999]

// 根部过渡点(占叶片长度比例)
root_transition = 0.33; // [0:0.01:1]

/* [分辨率] */
// 轮廓点数(0=自动)
profilepoints = 0; // [0:1:100]

// 叶片切片数(0=自动)
slices = 0; // [0:1:200]

/* [锥头选项] */
// 锥头类型
spinner_type = "parabolic"; // [none:无, parabolic:抛物线型, ogive:弧形]

// 锥头长度(毫米)
spinner_length = 20; // [5:1:50]

// 锥头直径(毫米)
spinner_diameter = 20; // [5:1:50]

// 弧形锥头鼻部半径(占直径比例)
ogive_noseradius = 0.2; // [0.05:0.01:0.24]

/* [隐藏参数] */
$fn = 100;

/* [显示选项] */
// 显示装配视图
show_assembly = true;

// 仅渲染单个叶片用于测试
show_single_blade = false;

// ...其余代码保持不变...
// Main assembly module
module mw_assembly_view() {
    if (show_assembly) {
        color("#99ccff") {
            if (show_single_blade) {
                propblade(propdia=propdia, hubdia=hubdia, bladepitch=bladepitch, 
                    maxchordfrac=maxchordfrac, hublen=hublen, elenfrac=elenfrac,
                    dir=dir, centerline=centerline, angle_sweep=angle_sweep,
                    te_thickness=te_thickness, fairing=fairing,
                    naca_root=naca_root, naca_mid=naca_mid, naca_tip=naca_tip,
                    root_transition=root_transition,
                    profilepoints=profilepoints, slices=slices);
            } else {
                difference() {
                    union() {
                        propeller(blades=blades, propdia=propdia, hubdia=hubdia, 
                            bladepitch=bladepitch, maxchordfrac=maxchordfrac,
                            hublen=hublen, elenfrac=elenfrac, dir=dir,
                            centerline=centerline, angle_sweep=angle_sweep,
                            te_thickness=te_thickness, fairing=fairing,
                            naca_root=naca_root, naca_mid=naca_mid, 
                            naca_tip=naca_tip, root_transition=root_transition,
                            profilepoints=profilepoints, slices=slices);
                        cylinder(h=hublen+2, d=hubdia);
                    }
                    translate([0,0,-1]) cylinder(hublen+2, d=holedia, $fn=36);
                }
            }
        }
        
        if (spinner_type != "none") {
            if (spinner_type == "parabolic") {
                paraboloc_spinner(spinner_length, spinner_diameter);
            } else if (spinner_type == "ogive") {
                ogive_spinner(spinner_length, spinner_diameter, ogive_noseradius);
            }
        }
    }
}

// Keep the original modules and functions below this line
// =======================================================================

// ---------- propeller ----------

// for argument explanation, see comments above or propblade() module below

module propeller(blades=2, propdia=100, hubdia=10, bladepitch=-1, maxchordfrac=0.15, hublen=10, elenfrac=1.1, dir=-1, centerline=0.3, angle_sweep=0, te_thickness=0.4, fairing=3, naca_root=9430, naca_mid=6412, naca_tip=4412, root_transition=0.33, profilepoints=0, slices=0) {
    for (a=[0:360/blades:359])
        rotate([0, 0, a])
            rotate([90, 0, 0])
                propblade(propdia=propdia, hubdia=hubdia, bladepitch=bladepitch, maxchordfrac=maxchordfrac, hublen=hublen, elenfrac=elenfrac, dir=dir, centerline=centerline, angle_sweep=angle_sweep, te_thickness=te_thickness, naca_root=naca_root, naca_mid=naca_mid, naca_tip=naca_tip, root_transition=root_transition, profilepoints=profilepoints, slices=slices, fairing=fairing);
}

// ---------- vertical blade ----------

module propblade(
    propdia=100,        // Propeller arc diameter.
    hubdia=10,          // Hub cylinder diameter where blade root starts.
    bladepitch=-1,      // Blade pitch - -1 defaults to 1/2 propdia.
    maxchordfrac=0.2,   // Maximum (at hub) chord length for zero attack if not wrapped around.
    hublen=10,          // Fore-to-aft length of hub cylinder.
    elenfrac=1.05,      // Fraction of blade length for ellipse profile.
    dir=-1,             // Rotation direction (1=CW, -1=CCW) when viewed from front.
    centerline=0.3,     // Align this fraction of each airfoil.
    angle_sweep=0,      // Blade sweep in degrees per millimeter.
    te_thickness = 0.4, // Thickness of airfoil trailing edge in mm.
    fairing=3,          // Millimeters to thicken (flair out) the root airfoil.
    naca_root = 8430,   // NACA code for airfoil at root (high-camber and thick).
    naca_mid = 6412,    // NACA code for main airfoil starting at transition point.
    naca_tip = 3412,    // NACA code for tip airfoil (less camber, more symmetric).
    root_transition = 0.33, // Fraction of blade length where root airfoil completes its transition to mid airfoil.
    profilepoints=0,    // Number of points for top and bottom of airfoils (0=calculate reasonable default).
    slices=0            // Slice blade length this many times (0=use 1mm slices).
    ) {

    // convert 4-digit NACA codes to dimensional values

    rootcamber = 0.01*floor(0.001*naca_root);       // root camber as fraction of chord
    rootcamberpos = 0.1*(floor(0.01*naca_root)%10); // root camber position as fraction of chord
    rootthickness = 0.01*(naca_root%100);           // root thickness as fraction of chord
    midcamber = 0.01*floor(0.001*naca_mid);         // mid camber as fraction of chord
    midcamberpos = 0.1*(floor(0.01*naca_mid)%10);   // mid camber position as fraction of chord
    midthickness = 0.01*(naca_mid%100);             // mid thickness as fraction of chord
    tipcamber = 0.01*floor(0.001*naca_tip);         // tip camber as fraction of chord
    tipcamberpos = 0.1*(floor(0.01*naca_tip)%10);   // tip camber position as fraction of chord
    tipthickness = 0.01*(naca_tip%100);             // tip thickness as fraction of chord

    // other initializations

    rhub = 0.5 * hubdia;            // hub radius
    length = 0.5*propdia - rhub;    // blade length from hub
    bpitch = bladepitch>=0 ? bladepitch : length;
    esemimajor = elenfrac*length;           // semimajor axis of ellipse
    maxchordlen = 2*maxchordfrac*length;    // full minor axis of ellipse
    hh = min(maxchordlen, hublen);  // hub height constraint for blade need not be longer than max chord
    ztrans = length*root_transition;    // transition point where root becomes mid airfoil    
    blen = length - ztrans;     // length of blade past the transition
    slicelen = slices > 0 ? length/slices : 1;  // size of cross-section slice
    airfoilsegs = profilepoints > 0 ?   // number of segments on one side of airfoil
        profilepoints : round(2*length * maxchordfrac);
    fairthickness =     // thickness adjustment for fairing at blade root
        fairing / ellipse_d(maxchordlen, hh, -dir*atan(bpitch/(2*PI*rhub)));

    // pxc is an array of airfoil profile cross-sections
    pxc = //fairing == 0 ? // make the blade if not making a fairing
        [
        // root section, from rhub to ztrans

        for(z=[(ztrans>0?0:-0.01*slicelen):slicelen:ztrans])
            let(
                rz = rhub + z,
                interp = sin(90*z/ztrans),
                fairthk = fairthickness * (fairing > 0 ? 1-sin(90*min(1,z/fairing)) : 0),
                thk = (midthickness-rootthickness)*interp+rootthickness + fairthk,
                cam = (midcamber-rootcamber)*interp+rootcamber,
                campos = (midcamberpos-rootcamberpos)*interp+rootcamberpos,
                attackangle = -dir*atan(bpitch/(2*PI*rz)),
                elen = maxchordlen*sqrt(1-z*z/(esemimajor*esemimajor)),
                chordlen = ellipse_d(elen, hh, attackangle))
                convert2dto3d(NACA_profile(airfoilsegs, cam, campos, thk, chordlen, centerline, dir, te_thickness), rz, attackangle, angle_sweep*sign(dir)),

        // main section, from ztrans to tip

        for(bz=[0:slicelen:blen+0.9*slicelen])
            let(z = min(bz,blen),
                ze = ztrans+z, rz = rhub + ze,
                interp = 1 - cos(90*z/blen),
                fairthk = fairthickness * (fairing > 0 ? 1-sin(90*min(1,ze/fairing)) : 0),
                thk = (tipthickness-midthickness)*interp+midthickness + fairthk,
                cam = (tipcamber-midcamber)*interp+midcamber,
                campos = (tipcamberpos-midcamberpos)*interp+midcamberpos,
                attackangle = -dir*atan(bpitch/(2*PI*rz)),
                elen = maxchordlen*sqrt(1-ze*ze/(esemimajor*esemimajor)),
                chordlen = ellipse_d(elen, hh, attackangle))
                convert2dto3d(NACA_profile(airfoilsegs, cam, campos, thk, chordlen, centerline, dir, te_thickness), rz, attackangle, angle_sweep*sign(dir))
        ];

    // assemble the blade object from all the airfoil cross-sections

    translate([0,(1-centerline)*hublen, 0])
        difference () {
            airfoil_polyhedron_stack(pxc); // connect all the profiles together into a big polyhedron
            // shave off anything extending under bottom of hub
            translate([0,-hublen*(1-centerline)-1,0]) cube([2*esemimajor+1, 2, 2*esemimajor+1], center=true);
        }

    // internal function to return chord length constrained by an ellipse defined by planform width and hub height
    function ellipse_d(majaxis, minaxis, angle) = // ellipse "diameter" at some angle
        let(a=0.5*majaxis, b=0.5*minaxis, bc = b*cos(angle), as = a*sin(angle))
            2 * a * b / sqrt(bc*bc + as*as);
}

// Given a stack of n-sided polygons in 3D space, connect them together into a polyhedron.
module airfoil_polyhedron_stack(stack) {
    nz = len(stack); // number of z layers
    np = len(stack[0]); // number of polygon vertices
    hnp = floor(np/2); // half of polygon vertices
    facets = [
        //[ for(j=[0:np-1]) j ], // close first opening
        for(k1=[0:hnp-1]) let(k2=k1+1, k3=np-2-k1, k4=np-1-k1) k1==k4 ? [k1,k2,k3] : [k1,k2,k3,k4],
        for(i=[0:nz-2])
            for(j=[0:np-1]) let(k1=i*np+j, k4=i*np+((j+1)%np), k2=k1+np, k3=k4+np)
                [k1, k2, k3, k4],
        // close last opening
        let(n=np*(nz-1)) for(k1=[0:hnp-1]) let(k2=k1+1, k3=np-2-k1, k4=np-1-k1) k2==k3 ? [k1+n,k4+n,k2+n] : [k1+n,k4+n,k3+n,k2+n]
    ];
    polyhedron(flatten(stack), facets, convexity=6);

    function flatten(v) = [ for (a=v) for (b=a) b ] ;
}

// ---------- spinners ----------

// parabolic propeller spinner
module paraboloc_spinner(length=20, diameter=20) {
    r = 0.5*diameter;
    p = [ [0,0], for(x=[0:0.05:1.001]) [ x*r, length*(1-x*x) ] ];
    rotate_extrude(angle=360, $fn=64) polygon(points=p);
}

// ogive (vertical slope base) propeller spinnner with rounded nose
// noseradius is a fraction of the diameter; must be <0.25
module ogive_spinner(length=20, diameter=20, noseradius=0.2) {
    rnose = noseradius*diameter;
    r = 0.5*diameter - rnose;
    ht = length-rnose;
    x = (ht*ht - r*r) / (2*r);
    circrad = x+r;
    astart = atan(ht/x);
    p = [ [0,rnose], for(a=[astart:-0.05*astart:-0.001]) [ circrad*cos(a)-x, circrad*sin(a) ] ];
    rotate_extrude(angle=360, $fn=64)
    difference() {
        offset(r=rnose, $fn=32) polygon(points=p);
        translate([-rnose-1,-1]) square(size=[rnose+1,length+2]);
        translate([-1,-rnose-1]) square(size=[r+2+rnose, rnose+1]);
    }
}

// ---------- demo modules ----------

// demo: three bladed propeller
module demo_nbladedprop(blades=3) {
    color("#99ccff") difference() {
        union() {
            propeller(blades=blades, propdia=100, hubdia=10, fairing=2);
            cylinder(11, d=10, $fn=60);
        }
        translate([0,0,-1]) cylinder(13, d=3, $fn=36);
    }
}

// demo: Random propeller with properties displayed in console.
module demo_random() {
    // set random values for blade parameters
    nblades = floor(rands(1,6,1)[0]) + 1;
    elenfrac = rands(1.01,1.8,1)[0];
    maxchordfrac = rands(0.15,0.9,1)[0];
    centerline = rands(0,1,1)[0]; // 1=trailing edge always on build plate
    sweep = rands(-1.5,1.5,1)[0];
    dia = 150;

    echo("Blades:", nblades);
    echo("Ellipse length fraction:", elenfrac);
    echo("Max chord fraction:", maxchordfrac);
    echo("Centerline at:", centerline, "of chord length");
    echo("Sweep angle degrees per radial millimeter:", sweep);

    // render the propeller
    propeller(blades=nblades, propdia=dia, hubdia=10, bladepitch=0.5*dia, hublen=20, elenfrac=elenfrac, maxchordfrac=maxchordfrac, centerline=centerline, angle_sweep=sweep);

    // put a spinner in the middle
    ogive_spinner(25,30);

    // randomly put an airfoil ring around the propeller
    if (rands(0,1,1)[0] < 0.4)
        rotate_extrude(angle=360, $fn=192) translate([75,0,0]) rotate([0,0,-90])
            polygon(NACA_profile(20, 0, 0.4, 0.15, chordlen=min(20,0.5*dia*maxchordfrac), origin=1, dir=1, te_thick=0.4));
}

// demo: 5 random props at once
module demo_5random() {
    for(a=[0:72:359]) rotate([0,0,a]) translate([140,0,0]) color("#99ccff") demo_random();
}

// demo: example boat propeller (blades only no hub)
module demo_boatpropblades() {
    hublen = 25;
    dia=80;
    translate([0,0,hublen+0.1]) rotate([180,0,0]) // flip upside-down
        propeller(blades=3, propdia=dia, hubdia=12, bladepitch=60, maxchordfrac=.7, hublen=hublen, elenfrac=1.001, naca_root=6410, naca_mid=4410, naca_tip=2410, centerline=0.5, angle_sweep=1);
}

// demo: printable boat propeller with spinner
module demo_printable_boatprop() {
    hublen = 25;
    dia=80;
    difference() {
        translate([0,0,hublen+0.1]) rotate([180,0,0]) // flip upside-down
            propeller(blades=3, propdia=dia, hubdia=6, bladepitch=60, maxchordfrac=.7, hublen=hublen, elenfrac=1.1, naca_root=6412, centerline=0, angle_sweep=2);
        translate([0,0,-6]) cylinder(6, d=dia+7);  // shave off leading edges that penetrate buildplate
    }
    cylinder(hublen-1, r=3, $fn=64);   // add the hub
    translate([0,0,hublen-1]) ogive_spinner(6, 6); // add the spinner
}

module demo_collection() {
    // 3-blade model airplane propeller
    translate([0,-30,0]) demo_nbladedprop(3);

    // 2-blde model airplane propeller
    translate([60,-35,0]) rotate([0,0,-70]) {
        propeller(blades=2, propdia=90, hubdia=8, bladepitch=-1, maxchordfrac=0.15, hublen=10, elenfrac=1.001, dir=-1, centerline=0.3, angle_sweep=0, profilepoints=32, slices=50);
        cylinder(11, r=4, $fn=32);
    }

    // toy boat propeller
    translate([0,30,0]) color("#99FF99") rotate([0,0,-8]) {
        demo_boatpropblades();
        cylinder(24, r=6, $fn=64);   // add the hub
        translate([0,0,24]) ogive_spinner(8, 12); // add the spinner
    }

    // toy rubber-band airplane propeller
    translate([-33,8,0]) color("#FF9999") rotate([0,0,60]) {
        propeller(blades=2, propdia=80, hubdia=3, bladepitch=70, maxchordfrac=0.5, hublen=6, elenfrac=1.001, dir=1, centerline=0.5, angle_sweep=0, profilepoints=32, slices=50);
        cylinder(7, r=1.5, $fn=30);
    }

    // computer cooling fan
    translate([62,18,0]) color("#999999") rotate([0,0,-5]) {
        propeller(blades=5, propdia=70, hubdia=30, bladepitch=70, maxchordfrac=0.75, hublen=16, elenfrac=2, dir=1, naca_root=6412, centerline=0.5, angle_sweep=0, profilepoints=32, slices=50, fairing=0);
        cylinder(16, r=15, $fn=30);
    }
}

// ---------- functions ----------

// convert 2D airfoil profile to 3D, tilting by attack angle and sweeping around the prop axle

function convert2dto3d(a, z, attackangle=0, anglesweep=0) = let(
    rd = 180/PI, n = len(a)-1,
    cosat = cos(attackangle), sinat = sin(attackangle),
    rotmatrix = [[cosat, -sinat], [sinat, cosat]],
    tilta = [ for(i=[0:n]) rotmatrix * a[i] ], // tilt the 2D airfoil
    sweeprot = z*anglesweep // angle position of airfoil origin
    ) [ // wrap the tilted 2D airfoil around the propeller shaft at distance z
        for(i=[0:n]) let(ang = rd*tilta[i][0]/z + sweeprot)
            [ z*sin(ang), tilta[i][1], z*cos(ang) ]
    ];

/*
From http://airfoiltools.com/airfoil/naca4digit
Given NACA 4-digit airfoil code ABCC
x = fraction of chord for which to calculate
M = maximum camber as fraction of chord (A/10)
P = position of maximum camber as fraction of chord (B/10)
T = thickness as a fraction of chord (CC/100)
*/
        
// primary function - return a polygon representing an airfoil cross-section
// chordlen = length of chord in mm
// origin = fraction of chord at which to center the polygon (0=front, 1=tail)
// dir = direction (1=facing left, -1=facing right)
// te_thick = trailing edge thickness
function NACA_profile(n, M, P, T, chordlen=1, origin=0, dir=1, te_thick=0.2) =
    dir < 0 ?
    [
        for (x=[1.0:-1/n:0.1/n]) NACA_profile_lower(x, M, P, T, chordlen, origin, dir, te_thick),
        for (x=[0:1/n:1-0.1/n]) NACA_profile_upper(x, M, P, T, chordlen, origin, dir, te_thick),
        NACA_profile_upper(1, M, P, T, chordlen, origin, dir, te_thick)
    ]
    : [
        for (x=[1.0:-1/n:0.1/n]) NACA_profile_upper(x, M, P, T, chordlen, origin, dir, te_thick),
        for (x=[0:1/n:1-0.1/n]) NACA_profile_lower(x, M, P, T, chordlen, origin, dir, te_thick),
        NACA_profile_lower(1, M, P, T, chordlen, origin, dir, te_thick)
    ];

// supporting functions

function NACA_camber(x, M, P) =
    let(a = 2*P*x-x*x, denom = x<P ? P : 1-P)
        M * (x<P ? a : 1-2*P+a) / (denom*denom);

function NACA_gradient(x, M, P) =
    let(num = 2*M*(P-x), denom = x<P ? P : 1-P)
        num / (denom*denom);

function NACA_thickness(x, T) =
    let(a0=0.2969, a1=-0.126, a2=-0.3516, a3=0.2843, a4=-0.1036 /*gives closed trailing edge; normally use -0.1015 to have a thickness*/ , x2=x*x)
        5*T*(a0*sqrt(x) + a1*x + a2*x2 + a3*x2*x + a4*x2*x2);

function NACA_profile_upper(x, M, P, T, chordlen=1, origin=0, dir=1, te_thick=0.2) = let(
    xp = 0.5 * (1 - cos(180*x)),
    theta = atan(NACA_gradient(xp, M, P)),
    yc = NACA_camber(xp, M, P),
    yt = NACA_thickness(xp, T))
        [chordlen*dir*(xp - yt * sin(theta) - origin), chordlen*(yc + yt * cos(theta)) + xp*te_thick/2];

function NACA_profile_lower(x, M, P, T, chordlen=1, origin=0, dir=1, te_thick=0.2) = let(
    xp = 0.5 * (1 - cos(180*x)),
    theta = atan(NACA_gradient(xp, M, P)),
    yc = NACA_camber(xp, M, P),
    yt = NACA_thickness(xp, T))
        [chordlen*dir*(xp + yt * sin(theta) - origin), chordlen*(yc - yt * cos(theta)) - xp*te_thick/2];

mw_assembly_view();

